home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
lin_alg.lha
/
lin_alg
/
matrix2.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
4KB
|
118 lines
// This may look like C code, but it is really -*- C++ -*-
/*
************************************************************************
*
* Linear Algebra Package
*
* Basic linear algebra operations, level 2
* Matrix transformations and multiplications
* of various types
*
************************************************************************
*/
#include "LinAlg.h"
#include <math.h>
#pragma implementation
/*
*------------------------------------------------------------------------
* Transpositions
*/
// Transpose a matrix
Matrix Matrix::transpose(void) const
return tm(ncols+col_lwb-1,nrows+row_lwb-1,col_lwb,row_lwb)
{
is_valid();
register REAL * rsp = elements; // Row source pointer
register REAL * tp = tm.elements;
// Matrix tm is traversed in the natural
// column-wise way, whilst the source matrix
// is scanned row-by-row
while( tp < tm.elements + tm.nelems )
{
register REAL * sp = rsp++; // sp = @ms[j,i] for i=0
while( sp < elements + nelems ) // Move tp to the next elem in the col,
*tp++ = *sp, sp += nrows; // sp to the next elem in the curr row
}
assert( tp == tm.elements + tm.nelems && rsp == elements + nrows );
}
/*
*------------------------------------------------------------------------
* General matrix multiplications
*/
// Compute target = target * source
// inplace,
// target being the matrix the operation is
// applied to.
Matrix& Matrix::operator *= (const Matrix& source)
{
is_valid();
source.is_valid();
if( row_lwb != source.col_lwb || ncols != source.nrows ||
col_lwb != source.col_lwb || ncols != source.ncols )
info(), source.info(),
_error("matrices above are unsuitable for the inplace multiplication");
REAL one_row[ncols]; // One row of the old_target matrix
register REAL * trp = elements; // Pointer to the i-th row
for(; trp < &elements[nrows]; trp++) // Go row-by-row in the target
{
register REAL *wrp, *orp; // work row pointers
for(wrp=trp,orp=&one_row[0]; orp < &one_row[ncols];)
*orp++ += *wrp, wrp += nrows; // Copy a row of old_target
register REAL *scp=source.elements; // source column pointer
for(wrp=trp; wrp < elements+nelems; wrp += nrows)
{
register double sum = 0; // Multiply a row of old_target
for(orp=one_row; orp < &one_row[ncols];) // by each col of source
sum += *orp++ * *scp++; // to get a row of new_target
*wrp = sum;
}
}
return *this;
}
// Create matrix C such that
// C = A * B
Matrix operator * (const Matrix& A, const Matrix& B)
return C(A.nrows,B.ncols,A.row_lwb,B.col_lwb)
{
A.is_valid();
B.is_valid();
if( A.ncols != B.nrows || A.row_lwb != B.col_lwb )
A.info(), B.info(),
_error("matrices above cannot be multiplied");
register REAL * arp; // Pointer to the i-th row of A
register REAL * bcp = B.elements; // Pointer to the j-th col of B
register REAL * cp = C.elements; // C is to be traversed in the natural
while( cp < C.elements + C.nelems ) // order, col-after-col
{
for(arp = A.elements; arp < A.elements + A.nrows; )
{
register double cij = 0;
while( arp < A.elements + A.nelems ) // Scan the i-th row of A and
cij += *bcp++ * *arp, arp += A.nrows; // the j-th col of B
*cp++ = cij;
bcp -= B.nrows; // Return bcp to the j-th col
arp -= A.nelems - 1; // arp points to (i+1)-th row
}
bcp += B.nrows; // We're done with j-th col of both
} // B and C. Set bcp to the (j+1)-th col
assert( cp == C.elements + C.nelems && bcp == B.elements + B.nelems );
}